This script visualizes and analyzes the mean amplitude of response
(meanbeta) to expected and unexpected events from three
“domains” - psychology-action, psychology-environment, and physics. Our
pre-registered confirmatory analyses center around the psychology-action
and physics videos.
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 1, digits = 3)
library(pacman)
pacman::p_load(tidyverse,
patchwork,
papaja,
ggforce,
conflicted,
here,
cowplot,
lme4,
lmerTest,
effects,
effectsize,
EMAtools,
performance,
parameters,
BayesFactor,
gt,
insight,
sjPlot,
lsmeans,
ggrepel,
GGally,
nptest,
boot,
boot.pval,
car)
conflict_prefer("select", "dplyr")
## [conflicted] Will prefer dplyr::select over any other package
conflict_prefer("filter", "dplyr")
## [conflicted] Will prefer dplyr::filter over any other package
conflict_prefer("rename", "dplyr")
## [conflicted] Will prefer dplyr::rename over any other package
conflict_prefer("mutate", "dplyr")
## [conflicted] Will prefer dplyr::mutate over any other package
conflict_prefer("summarise", "dplyr")
## [conflicted] Will prefer dplyr::summarise over any other package
conflict_prefer("lmer", "lmerTest")
## [conflicted] Will prefer lmerTest::lmer over any other package
conflict_prefer("here", "here")
## [conflicted] Will prefer here::here over any other package
source(here("helper_funs.R"))
options(contrasts = c("contr.sum", "contr.poly"))
theme_set(theme_cowplot(10))
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] car_3.0-8 boot.pval_0.4 boot_1.3-25
## [4] nptest_1.0-3 GGally_2.1.2 ggrepel_0.8.2
## [7] lsmeans_2.30-0 emmeans_1.4.8 sjPlot_2.8.10
## [10] insight_0.15.0 gt_0.8.0 BayesFactor_0.9.12-4.4
## [13] coda_0.19-3 parameters_0.16.0 performance_0.8.0
## [16] EMAtools_0.1.3 effectsize_0.6.0.1 effects_4.1-4
## [19] carData_3.0-4 lmerTest_3.1-2 lme4_1.1-23
## [22] Matrix_1.2-18 cowplot_1.0.0 here_1.0.1
## [25] conflicted_1.1.0 ggforce_0.3.4 papaja_0.1.0.9997
## [28] tinylabels_0.2.3 patchwork_1.1.2 forcats_0.5.0
## [31] stringr_1.5.0 dplyr_1.0.9 purrr_0.3.4
## [34] readr_1.3.1 tidyr_1.2.0 tibble_3.1.7
## [37] ggplot2_3.4.1 tidyverse_1.3.0 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.4.1 plyr_1.8.6
## [4] splines_4.0.2 TH.data_1.0-10 digest_0.6.31
## [7] htmltools_0.5.4 fansi_1.0.3 magrittr_2.0.3
## [10] memoise_2.0.1 openxlsx_4.1.5 modelr_0.1.8
## [13] sandwich_2.5-1 colorspace_2.0-3 blob_1.2.1
## [16] rvest_0.3.6 mitools_2.4 rbibutils_2.2.8
## [19] haven_2.4.3 xfun_0.36 crayon_1.5.1
## [22] jsonlite_1.8.4 survival_3.1-12 zoo_1.8-8
## [25] glue_1.6.2 polyclip_1.10-0 gtable_0.3.0
## [28] MatrixModels_0.4-1 sjstats_0.18.0 sjmisc_2.8.5
## [31] abind_1.4-5 scales_1.2.0 mvtnorm_1.1-1
## [34] DBI_1.1.2 ggeffects_0.16.0 Rcpp_1.0.8
## [37] xtable_1.8-4 foreign_0.8-80 survey_4.0
## [40] datawizard_0.2.3 httr_1.4.2 RColorBrewer_1.1-3
## [43] ellipsis_0.3.2 pkgconfig_2.0.3 reshape_0.8.8
## [46] farver_2.1.0 nnet_7.3-14 sass_0.4.4
## [49] dbplyr_1.4.4 utf8_1.2.2 tidyselect_1.1.2
## [52] rlang_1.0.6 munsell_0.5.0 cellranger_1.1.0
## [55] tools_4.0.2 cachem_1.0.6 cli_3.5.0
## [58] generics_0.1.2 sjlabelled_1.1.7 broom_0.8.0
## [61] evaluate_0.19 fastmap_1.1.0 yaml_2.3.6
## [64] knitr_1.41 fs_1.5.2 zip_2.2.0
## [67] pbapply_1.6-0 nlme_3.1-148 xml2_1.3.2
## [70] compiler_4.0.2 rstudioapi_0.13 curl_4.3
## [73] reprex_0.3.0 statmod_1.4.34 tweenr_1.0.2
## [76] bslib_0.4.2 stringi_1.7.8 lattice_0.20-41
## [79] nloptr_1.2.2.2 vctrs_0.5.1 pillar_1.7.0
## [82] lifecycle_1.0.3 Rdpack_2.3.1 jquerylib_0.1.4
## [85] estimability_1.3 data.table_1.13.0 R6_2.5.1
## [88] rio_0.5.16 codetools_0.2-16 MASS_7.3-51.6
## [91] assertthat_0.2.1 rprojroot_2.0.3 withr_2.5.0
## [94] multcomp_1.4-13 bayestestR_0.11.5 hms_0.5.3
## [97] DataCombine_0.2.21 grid_4.0.2 minqa_1.2.4
## [100] rmarkdown_2.19.2 numDeriv_2016.8-1.1 lubridate_1.7.9
regioninfo <- read.csv(here("input_data/manyregions_info.csv"))
DT::datatable(regioninfo %>% dplyr::arrange(ROI_category, focal_region), options = list(scrollX = TRUE))
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
univariate_data <-
read_rds(here("outputs/univariate_data/study2_univariate_data_allregions_alltopN.Rds"))
univariate_data$event <- relevel(univariate_data$event, ref = "unexp") # for some reason this is required for the model coefficients to reflect the intuitive def that unexpected > expected is a positive coefficient
univariate_summary_domain <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain.Rds"))
univariate_summary_domain_runs12 <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12.Rds"))
univariate_summary_domain_allruns <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_allruns.Rds"))
focal_data <- univariate_data %>%
filter(focal_region == 1)
focal_summary_domain <- univariate_summary_domain %>%
filter(focal_region == 1)
focal_summary_domain_runs12 <- univariate_summary_domain_runs12 %>%
filter(focal_region == 1)
# focal_summary_task_runs12 <- univariate_summary_task_runs12 %>%
# filter(focal_region == 1)
regions <- levels(as.factor(focal_data$ROI_name_final))
First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).
plot_univar_runbyrun("superiortemporal_L", -1, 4) | plot_univar_runbyrun("superiortemporal_R", -1, 4)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_L", -1, 4) | plot_univar_runbyrun("supramarginal_R", -1, 4)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("antParietal_bilateral", -1, 6) | plot_univar_runbyrun("precentral_inferiorfrontal_R", -1, 6)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("V1_bilateral", -1, 11) | plot_univar_runbyrun("MT_bilateral", -1, 11)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
These plots go into Figure 3 of the main text.
# DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
# plot_univar_runs12("supramarginal", 0, 4)
#
# DG_runs12_plots <- plot_univar_runs12("MD", 0, 6) | plot_univar_runs12("early_visual", 0, 11)
#
# DS_runs12_plots / DG_runs12_plots
(plot_univar_runs12("superiortemporal_L", 0, 4) | plot_univar_runs12("superiortemporal_R", 0, 4) )
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
(plot_univar_runs12("supramarginal_L", 0, 4) | plot_univar_runs12("supramarginal_R", 0, 4))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
(plot_univar_runs12("antParietal_bilateral", 0, 6) | plot_univar_runs12("precentral_inferiorfrontal_R", 0, 6))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
summary_domain_allruns <- read_rds(here("outputs/descriptive_summaries/study2_univariate_summary_domain_allruns.Rds"))
Univariate response to psych_env events, over all runs, for all focal regions.
STS <- plot_both_allruns("superiortemporal_", -1, 10)
SMG <- plot_both_allruns("supramarginal_", -1, 10)
V1 <- plot_both_allruns("V1_bilateral", -1, 10)
MT <- plot_both_allruns("MT_bilateral", -1, 10)
APC <- plot_both_allruns("antParietal_bilateral", -1, 10)
RFC <- plot_both_allruns("precentral_inferiorfrontal_R", -1, 10)
STS | SMG | V1 | MT | APC | RFC + plot_layout(widths = c(2, 2, 1, 1, 1, 1))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
study2_univariate_summary_task_singlebetas <- readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_task_singlebetas.Rds")) %>%
filter(focal_region == "1")
regions <- levels(as.factor(study2_univariate_summary_task_singlebetas$ROI_name_final))
plot_univar <- function(region, ymin, ymax) {
plotobject <- ggplot(data = study2_univariate_summary_task_singlebetas %>%
filter(str_detect(ROI_name_final, region)),
aes(x = event, y = meanbeta, fill = task)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint")) + extracted_run_number, nrow = 2) +
# scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}
plot_univar("superiortemporal_L", -2, 3)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("superiortemporal_R", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("supramarginal_L", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("supramarginal_R", -2, 4)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("V1_bilateral", -1, 11)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("MT_bilateral", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("antParietal_bilateral", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
plot_univar("precentral_inferiorfrontal_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
Following the preregistration plan, we check whether we are able to take the second level cope values (averaging across 4 runs) as the main dependent measure, or whether we need to restrict the analysis to just the first 2 runs.
First we fit a model including an interaction between run number and event.
focal_data_checking <- focal_data %>%
filter(event != "fam",
top_n_voxels == "100",
domain != "both",
str_length(extracted_run_number) == 4) %>%
filter(ROI_name_final != "V1_bilateral",
ROI_name_final != "MT_bilateral")
focal_data_checking$fixation_position <- as.factor(focal_data_checking$fixation_position)
focal_data_checking$event <- relevel(focal_data_checking$event, ref = "unexp")
checkmodel <- lmer(
data = focal_data_checking,
formula = meanbeta ~ extracted_run_number * event + (1|subjectID) + (1|ROI_name_final)
)
checkmodel %>%
apa_print() %>%
apa_table()
| Term | \(\hat{\beta}\) | 95% CI | \(t\) | \(\mathit{df}\) | \(p\) |
|---|---|---|---|---|---|
| Intercept | 2.20 | [1.25, 3.16] | 4.51 | 6.96 | .003 |
| Extracted run number1 | 0.41 | [0.27, 0.55] | 5.83 | 3,004.27 | < .001 |
| Extracted run number2 | 0.15 | [0.02, 0.29] | 2.18 | 3,004.27 | .030 |
| Extracted run number3 | -0.15 | [-0.29, -0.01] | -2.06 | 3,006.55 | .039 |
| Event1 | 0.09 | [0.01, 0.17] | 2.27 | 3,003.97 | .023 |
| Extracted run number1 \(\times\) Event1 | 0.06 | [-0.08, 0.19] | 0.78 | 3,003.97 | .433 |
| Extracted run number2 \(\times\) Event1 | 0.07 | [-0.07, 0.21] | 1.00 | 3,003.97 | .318 |
| Extracted run number3 \(\times\) Event1 | -0.07 | [-0.21, 0.07] | -0.92 | 3,003.97 | .356 |
summary(checkmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: focal_data_checking
##
## REML criterion at convergence: 13783
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.973 -0.603 -0.053 0.550 5.440
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.19 1.09
## ROI_name_final (Intercept) 1.20 1.10
## Residual 5.14 2.27
## Number of obs: 3048, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2045 0.4890 6.9622 4.51 0.0028
## extracted_run_number1 0.4138 0.0710 3004.2738 5.83 6.1e-09
## extracted_run_number2 0.1544 0.0710 3004.2738 2.18 0.0297
## extracted_run_number3 -0.1484 0.0720 3006.5502 -2.06 0.0393
## event1 0.0934 0.0411 3003.9740 2.27 0.0230
## extracted_run_number1:event1 0.0556 0.0709 3003.9740 0.78 0.4332
## extracted_run_number2:event1 0.0709 0.0709 3003.9740 1.00 0.3180
## extracted_run_number3:event1 -0.0662 0.0717 3003.9740 -0.92 0.3557
##
## (Intercept) **
## extracted_run_number1 ***
## extracted_run_number2 *
## extracted_run_number3 *
## event1 *
## extracted_run_number1:event1
## extracted_run_number2:event1
## extracted_run_number3:event1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1 -0.001
## extrctd_r_2 -0.001 -0.329
## extrctd_r_3 0.002 -0.338 -0.338
## event1 0.000 0.000 0.000 0.000
## extrct__1:1 0.000 0.000 0.000 0.000 -0.005
## extrct__2:1 0.000 0.000 0.000 0.000 -0.005 -0.330
## extrct__3:1 0.000 0.000 0.000 0.000 0.014 -0.337 -0.337
lsmeans(checkmodel, pairwise ~ event | extracted_run_number)$contrasts
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.
## extracted_run_number = run1:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.298 0.164 Inf 1.821 0.0690
##
## extracted_run_number = run2:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.328 0.164 Inf 2.008 0.0450
##
## extracted_run_number = run3:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.054 0.166 Inf 0.327 0.7440
##
## extracted_run_number = run4:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.066 0.164 Inf 0.405 0.6850
##
## Degrees-of-freedom method: asymptotic
model_habituation_results <- cbind(gen.m(checkmodel), gen.ci(checkmodel)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(checkmodel, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()
We then compare this model to a simpler model without this interaction.
checkmodel_simpler <- lmer(
data = focal_data_checking,
formula = meanbeta ~ event + extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)
summary(checkmodel_simpler)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event + extracted_run_number + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: focal_data_checking
##
## REML criterion at convergence: 13774
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.947 -0.602 -0.057 0.553 5.414
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.19 1.09
## ROI_name_final (Intercept) 1.20 1.10
## Residual 5.14 2.27
## Number of obs: 3048, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2045 0.4890 6.9620 4.51 0.0028 **
## event1 0.0939 0.0411 3006.9741 2.29 0.0222 *
## extracted_run_number1 0.4138 0.0710 3007.2741 5.83 6.1e-09 ***
## extracted_run_number2 0.1544 0.0710 3007.2741 2.18 0.0297 *
## extracted_run_number3 -0.1484 0.0720 3009.5523 -2.06 0.0393 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 ext__1 ext__2
## event1 0.000
## extrctd_r_1 -0.001 0.000
## extrctd_r_2 -0.001 0.000 -0.329
## extrctd_r_3 0.002 0.000 -0.338 -0.338
# checkmodel_simpler %>%
# apa_print() %>%
# apa_table()
anova(checkmodel, checkmodel_simpler)
## refitting model(s) with ML (instead of REML)
knitr::kable(model_habituation_results)
| B | SE | df | t | p | CI_95_lower | CI_95_upper | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 2.204 | 0.489 | 6.96 | 4.508 | 0.003 | 2.208 | 2.323 |
| extracted_run_number1 | 0.414 | 0.071 | 3004.27 | 5.831 | 0.000 | 1.183 | 3.226 |
| extracted_run_number2 | 0.154 | 0.071 | 3004.27 | 2.175 | 0.030 | 0.275 | 0.553 |
| extracted_run_number3 | -0.148 | 0.072 | 3006.55 | -2.062 | 0.039 | 0.015 | 0.293 |
| event1 | 0.093 | 0.041 | 3003.97 | 2.274 | 0.023 | -0.289 | -0.007 |
| extracted_run_number1:event1 | 0.056 | 0.071 | 3003.97 | 0.784 | 0.433 | 0.013 | 0.174 |
| extracted_run_number2:event1 | 0.071 | 0.071 | 3003.97 | 0.999 | 0.318 | -0.083 | 0.195 |
| extracted_run_number3:event1 | -0.066 | 0.072 | 3003.97 | -0.924 | 0.356 | -0.068 | 0.210 |
knitr::kable(model_habituation_results_byrun)
| contrast | extracted_run_number | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| unexp - exp | run1 | 0.298 | 0.164 | 3004 | 1.821 | 0.069 |
| unexp - exp | run2 | 0.328 | 0.164 | 3004 | 2.008 | 0.045 |
| unexp - exp | run3 | 0.054 | 0.166 | 3004 | 0.327 | 0.744 |
| unexp - exp | run4 | 0.066 | 0.164 | 3004 | 0.405 | 0.685 |
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))
This is a grey area. On the one hand we do not observe a significant run x event interaction. On the other hand, the primary VOE effect (unexp > expected) is indeed variable across runs, and stronger in runs 1 and 2 than runs 3 and 4, when considering psychology and physics events). So on balance, we have the grounds for focusing the analysis on runs 1-2.
Here we fit a model including main effects of domain and event (domain-general and domain-specific regions) and also their interaction (domain-specific regions only), and organize and print the outputs.
NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean. This also includes standardized betas (effect sizes in SDs.)
data_for_modeling <- univariate_data %>%
filter(event != "fam",
# top_n_voxels == "100",
!str_detect(extracted_run_number, "fold"),
!str_detect(extracted_run_number, "all_runs"),
domain != "both")
ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
top_n_voxels <- unique(data_for_modeling$top_n_voxels)
# top_n_voxels <- 100
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "event:domain")
modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun = run
if (str_length(curRun) == 4) { # run1, run2, etc
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
} else if (curRun == "all_runs") {
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
str_length(extracted_run_number) == 4)
} else if (curRun == "runs12") {
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN) %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")
}
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
} else {
ROI_category <- "general"
}
# had to remove run random int (convergence))
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID))
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}
}
}
modelsummaries <- modelsummaries_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries$effect <-
as.factor(modelsummaries$effect)
levels(modelsummaries$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(modelsummaries) <- NULL
# write_rds(modelsummaries, here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))
modelsummaries <- read_rds(here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))
DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))
Below we check, for each domain-specific region (left and right STS and SMG), whether a model including an interaction between domain and event fits better than a model without the interaction.
focal_data_100_runs12 <- focal_data %>%
filter(top_n_voxels == "100",
domain != "both",
event != "fam") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")
LSMG_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + (1|subjectID)
)
plot(allEffects(LSMG_model))
# check_model(LSMG_model)
check_outliers(LSMG_model)
## OK: No outliers detected.
summary(LSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | subjectID)
## Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L")
##
## REML criterion at convergence: 948
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4287 -0.5510 0.0345 0.5874 2.6185
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.02 2.00
## Residual 1.58 1.26
## Number of obs: 256, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3899 0.3630 31.0000 6.58 2.4e-07 ***
## event1 0.0672 0.0787 221.0000 0.85 0.3942
## domain1 0.2580 0.0787 221.0000 3.28 0.0012 **
## event1:domain1 0.1877 0.0787 221.0000 2.39 0.0178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000
LSMG_results_bydomain_top100 <-
lsmeans(LSMG_model, pairwise ~ event |
domain)$contrasts %>% as.data.frame()
knitr::kable(LSMG_results_bydomain_top100)
| contrast | domain | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| unexp - exp | physics | 0.510 | 0.222 | 221 | 2.29 | 0.023 |
| unexp - exp | psychology | -0.241 | 0.222 | 221 | -1.08 | 0.280 |
# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))
LSMG_model_simpler <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event + domain + (1|subjectID)
)
# check_model(LSMG_model_simpler)
check_outliers(LSMG_model_simpler)
## OK: No outliers detected.
LSMG_model_comparison <- anova(LSMG_model_simpler, LSMG_model) %>%
as.data.frame()
## refitting model(s) with ML (instead of REML)
# write_rds(LSMG_model_comparison, path = here("outputs/univariate_results/LSMG_model_comparison_top100.Rds"))
RSMG_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(RSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R")
##
## REML criterion at convergence: 996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.437 -0.585 0.036 0.515 3.808
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 6.484 2.546
## extracted_run_number (Intercept) 0.098 0.313
## Residual 1.829 1.353
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.5420 0.5087 14.6472 5.00 0.00017 ***
## event1 0.1613 0.0845 219.9999 1.91 0.05767 .
## domain1 0.5037 0.0845 219.9999 5.96 1e-08 ***
## event1:domain1 0.0566 0.0845 219.9999 0.67 0.50410
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000
#
# RSMG_model %>%
# apa_print() %>%
# apa_table()
# check_model(RSMG_model)
check_outliers(RSMG_model)
## OK: No outliers detected.
lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.436 0.239 220 1.822 0.0700
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.209 0.239 220 0.876 0.3820
##
## Degrees-of-freedom method: kenward-roger
RSTS_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(RSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data:
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R")
##
## REML criterion at convergence: 960
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4296 -0.5228 0.0122 0.5633 3.0469
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.6191 1.618
## extracted_run_number (Intercept) 0.0135 0.116
## Residual 1.7596 1.326
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0199 0.3090 20.7114 6.54 1.9e-06 ***
## event1 0.1234 0.0829 220.0000 1.49 0.1381
## domain1 -0.2421 0.0829 220.0000 -2.92 0.0039 **
## event1:domain1 0.0539 0.0829 220.0000 0.65 0.5160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000
# RSTS_model %>%
# apa_print() %>%
# apa_table()
# check_model(RSTS_model)
check_outliers(RSTS_model)
## OK: No outliers detected.
lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.355 0.234 220 1.512 0.1320
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.139 0.234 220 0.592 0.5540
##
## Degrees-of-freedom method: kenward-roger
LSTS_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(LSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data:
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L")
##
## REML criterion at convergence: 1041
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.539 -0.593 0.035 0.503 3.454
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.337 2.083
## extracted_run_number (Intercept) 0.142 0.376
## Residual 2.362 1.537
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.6181 0.4644 6.5575 1.33 0.228
## event1 0.0624 0.0960 220.0000 0.65 0.517
## domain1 -0.2998 0.0960 220.0000 -3.12 0.002 **
## event1:domain1 -0.0658 0.0960 220.0000 -0.69 0.494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000
# LSTS_model %>%
# apa_print() %>%
# apa_table()
# check_model(LSTS_model)
check_outliers(LSTS_model)
## Warning: 1 outliers detected (cases 52).
outliers_list <- check_outliers(LSTS_model)
insight::get_data(LSTS_model)[outliers_list, ]
LSTS_model_nooutliers <- lmer(
data = insight::get_data(LSTS_model) %>%
filter(meanbeta != -11.6),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
# summary(LSTS_model_nooutliers)
LSTS_model_nooutliers %>%
apa_print() %>%
apa_table()
| Term | \(\hat{\beta}\) | 95% CI | \(t\) | \(\mathit{df}\) | \(p\) |
|---|---|---|---|---|---|
| Intercept | 0.62 | [-0.29, 1.53] | 1.33 | 6.56 | .228 |
| Event1 | 0.06 | [-0.13, 0.25] | 0.65 | 220.00 | .517 |
| Domain1 | -0.30 | [-0.49, -0.11] | -3.12 | 220.00 | .002 |
| Event1 \(\times\) Domain1 | -0.07 | [-0.25, 0.12] | -0.69 | 220.00 | .494 |
lsmeans(LSTS_model_nooutliers, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.0069 0.272 220 -0.026 0.9800
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.2564 0.272 220 0.944 0.3460
##
## Degrees-of-freedom method: kenward-roger
focal_modelsummaries <- modelsummaries %>%
filter(focal_region == 1)
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs"))
ggplot(focal_modelsummaries %>%
filter(str_length(extracted_run_number) == 4,
effect != "intercept"), aes(extracted_run_number, B_standardized, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Standardized Beta") +
facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral")), nrow = 3) +
# coord_cartesian(ylim = c(-.75, .75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Then we repeat the same modeling procedure on these regions. But this time, we fit an interaction between event and domain in all regions, because some MD regions appear to have a preference for physical events, and it is an open question, if this is true, whether they too could encode physical prediction error.
DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect))
This set of regions includes left and right V1, MT, insula, and IFG (rather than bilateral). The bilateral versions of these ROIs are not included in this portion of the analysis.
DT::datatable(modelsummaries %>% filter(focal_region == 0, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))
NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean. This includes effect sizes (e.g. 0.2 = 0.2 standard deviations from the grand mean; for a two-level factor, 0.4 standard deviations from the reference level)
focal_data_both_100 <- focal_data %>%
filter(event != "fam",
top_n_voxels == "100",
domain == "both")
options(contrasts = c("contr.sum", "contr.poly"))
checkdata_both <- focal_data_both_100 %>%
filter(ROI_name_final != "V1_bilateral",
ROI_name_final != "MT_bilateral",
str_length(extracted_run_number) == 4)
checkdata_both$fixation_position <- as.factor(checkdata_both$fixation_position)
checkdata_both$event <- relevel(checkdata_both$event, ref = "unexp")
checkmodel_bothdomains <- lmer(
data = checkdata_both,
formula = meanbeta ~ event * extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)
summary(checkmodel_bothdomains)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * extracted_run_number + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: checkdata_both
##
## REML criterion at convergence: 7354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.428 -0.547 -0.028 0.538 5.262
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.11 1.05
## ROI_name_final (Intercept) 1.34 1.16
## Residual 6.83 2.61
## Number of obs: 1524, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0973 0.5121 6.6088 4.10 0.0052
## event1 0.2175 0.0669 1479.8794 3.25 0.0012
## extracted_run_number1 0.5830 0.1157 1480.2519 5.04 5.2e-07
## extracted_run_number2 0.1838 0.1157 1480.2519 1.59 0.1123
## extracted_run_number3 -0.4030 0.1172 1483.0707 -3.44 0.0006
## event1:extracted_run_number1 -0.1024 0.1156 1479.8794 -0.89 0.3760
## event1:extracted_run_number2 -0.0745 0.1156 1479.8794 -0.64 0.5194
## event1:extracted_run_number3 0.2402 0.1169 1479.8794 2.06 0.0400
##
## (Intercept) **
## event1 **
## extracted_run_number1 ***
## extracted_run_number2
## extracted_run_number3 ***
## event1:extracted_run_number1
## event1:extracted_run_number2
## event1:extracted_run_number3 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 ext__1 ext__2 ext__3 e1:__1 e1:__2
## event1 0.000
## extrctd_r_1 -0.001 0.000
## extrctd_r_2 -0.001 0.000 -0.329
## extrctd_r_3 0.002 0.000 -0.338 -0.338
## evnt1:xt__1 0.000 -0.005 0.000 0.000 0.000
## evnt1:xt__2 0.000 -0.005 0.000 0.000 0.000 -0.330
## evnt1:xt__3 0.000 0.014 0.000 0.000 0.000 -0.337 -0.337
# checkmodel %>%
# apa_print() %>%
# apa_table()
plot(allEffects(checkmodel_bothdomains))
lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number)$contrasts
## extracted_run_number = run1:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.230 0.267 1480 0.860 0.3880
##
## extracted_run_number = run2:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.286 0.267 1480 1.070 0.2840
##
## extracted_run_number = run3:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.915 0.271 1480 3.380 0.0010
##
## extracted_run_number = run4:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.308 0.267 1480 1.160 0.2480
##
## Degrees-of-freedom method: kenward-roger
model_habituation_results <- cbind(gen.m(checkmodel_bothdomains), gen.ci(checkmodel_bothdomains)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_psych-env_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_psych-env_perrun.Rds"))
univariate_data_both <- univariate_data %>%
filter(event != "fam",
# top_n_voxels == "100",
domain == "both",
str_length(extracted_run_number) == 4)
ROIs <- levels(as.factor(univariate_data_both$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
# top_n_voxels <- 100
modelsummaries_both_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun <- run
if (curRun != "all_runs" & curRun != "runs12") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
} else if (curRun == "all_runs") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN)
} else if (curRun == "runs12") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
}
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
}
else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID))
focal_region <- ROIdata$focal_region[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4,],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
modelsummaries_both_init <-
rbind(modelsummaries_both_init, modelsummary)
}
}
}
modelsummaries_both <- modelsummaries_both_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries_both$effect <-
as.factor(modelsummaries_both$effect)
levels(modelsummaries_both$effect) <-
c("intercept", "event")
rownames(modelsummaries_both) <- NULL
# write_rds(modelsummaries_both, here("outputs/univariate_results/univar_results_psych-env_allregions_alltopN.Rds"))
modelsummaries_both <- read_rds( here("outputs/univariate_results/univar_results_psych-env_allregions_alltopN.Rds"))
Effects by run and ROI size
modelsummaries_both$extracted_run_number <- factor(modelsummaries_both$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs"))
modelsummaries_both_focal <- modelsummaries_both %>% filter(focal_region == 1, effect == "event")
ggplot(modelsummaries_both_focal %>% filter(str_length(extracted_run_number) == 4), aes(extracted_run_number, B_standardized, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
# geom_errorbar(aes(ymin = B - SE, ymax = B + SE, width = 0.1)) +
# ylab("Bayes Factor (log)") +
ylab("Standardized Beta") +
facet_wrap(~ effect + ROI, nrow = 1) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
focal_modelsummaries <- modelsummaries_both %>%
filter(focal_region == 1,
str_length(extracted_run_number) == 4)
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4"))
ggplot(
focal_modelsummaries %>%
filter(str_length(extracted_run_number) == 4,
effect == "event"),
aes(extracted_run_number, B_standardized, colour = top_n_voxels)
) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Standardized Beta") +
facet_wrap( ~ effect + factor(
ROI,
levels = c(
"supramarginal_L",
"supramarginal_R",
"superiortemporal_L",
"superiortemporal_R",
"antParietal_bilateral",
"precentral_inferiorfrontal_R",
"V1_bilateral",
"MT_bilateral"
)
), nrow = 1) +
# coord_cartesian(ylim = c(-.75, .75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))
DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "all_runs") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))
DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "runs12") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))
Here we repeat the same analysis as confirmatory analyses above, except that we use VOE runs 3-4 to localize STS and SMG, and we search for voxels in the domain general ROI parcels that respond more to unexp than exp, to ask whether there are any voxels that show the VOE response in those regions, without considering the localizer data. NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean. This includes effect sizes (e.g. 0.2 = 0.2 standard deviations from the grand mean; for a two-level factor, 0.4 standard deviations from the reference level)
options(contrasts = c("contr.sum", "contr.poly"))
univariate_data_VOEextract <-
readRDS(here("outputs", "univariate_data", "study2_univariate_data_VOEextract.Rds"))
univariate_summary_domain_VOEextract <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_VOEextract.Rds"))
univariate_summary_domain_runs12_VOEextract <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12_VOEextract.Rds"))
focal_data_VOEextract <- univariate_data_VOEextract %>%
filter(focal_region == 1)
focal_summary_domain_VOEextract <- univariate_summary_domain_VOEextract %>%
filter(focal_region == 1)
focal_summary_domain_runs12_VOEextract <- univariate_summary_domain_runs12_VOEextract %>%
filter(focal_region == 1)
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
regions <- levels(as.factor(focal_data_VOEextract$ROI_name_final))
First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).
plot_univar_runbyrun("superiortemporal_L")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("superiortemporal_R")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_L")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_R")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("antParietal_bilateral")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("precentral_inferiorfrontal_R")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("V1_bilateral")
## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("MT_bilateral")
## Warning: Using alpha for a discrete variable is not advised.
DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
plot_univar_runs12("supramarginal", 0, 4)
DG_runs12_plots <- plot_univar_runs12("precentral_inferiorfrontal_R", 0, 4) | plot_univar_runs12("antParietal_bilateral", 0, 4) | plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)
DS_runs12_plots / DG_runs12_plots
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
focal_data_100_runs12_VOEextract <- focal_data_VOEextract %>%
filter(
event != "fam",
top_n_voxels == "100",
domain != "both",
(extracted_run_number == "run1" | extracted_run_number == "run2"))
ROIs <- levels(as.factor(focal_data_100_runs12_VOEextract$ROI_name_final))
focal_data_100_runs12_VOEextract$event <- relevel(focal_data_100_runs12_VOEextract$event, ref = "exp")
focal_modelsummaries_100_VOEextract_init <- data.frame()
options(contrasts = c("contr.sum", "contr.poly"))
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
focal_data_100_runs12_VOEextract %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID)) # had to remove run random int (convergence))
} else {
ROI_category <- "general"
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID))
}
model_just_maineffects <- update(model, formula = ~ . -event:domain) # Without interaction term
model_just_domain <- update(model_just_maineffects, formula = ~ . -event)
model_just_event <- update(model_just_maineffects, formula = ~ . -domain)
BF_BIC_interaction <- data.frame(BF = exp((BIC(model_just_maineffects) - BIC(model))/2), BF_interpretation = interpret_bf(exp((BIC(model_just_maineffects) - BIC(model))/2))) # BICs to Bayes factor
BF_BIC_event <- data.frame(BF = exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2),
BF_interpretation = interpret_bf(exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2)))
BF_BIC_domain <- data.frame(BF = exp((BIC(model_just_event) - BIC(model_just_maineffects))/2),
BF_interpretation = interpret_bf(exp((BIC(model_just_event) - BIC(model_just_maineffects))/2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
focal_modelsummaries_100_VOEextract_init <-
rbind(focal_modelsummaries_100_VOEextract_init, modelsummary)
}
focal_modelsummaries_100_VOEextract <- focal_modelsummaries_100_VOEextract_init %>%
mutate_at(c(5:13), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
focal_modelsummaries_100_VOEextract$effect <-
as.factor(focal_modelsummaries_100_VOEextract$effect)
levels(focal_modelsummaries_100_VOEextract$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(focal_modelsummaries_100_VOEextract) <- NULL
# write_rds(focal_modelsummaries_100_VOEextract, here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))
focal_modelsummaries_100_VOEextract <- read_rds(here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))
DT::datatable(dplyr::arrange(focal_modelsummaries_100_VOEextract, effect), options = list(scrollX = TRUE))
focal_data_100_runs1234_VOEextract_bothdomains <-
focal_data_VOEextract %>%
filter(
event != "fam",
top_n_voxels == "100",
domain == "both",
str_length(extracted_run_number) == 4
)
ROIs <-
levels(as.factor(
focal_data_100_runs1234_VOEextract_bothdomains$ROI_name_final
))
focal_data_100_runs1234_VOEextract_bothdomains$event <-
relevel(focal_data_100_runs1234_VOEextract_bothdomains$event, ref = "unexp")
focal_modelsummaries_100_VOEextract_bothdomains_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
cur_top_voxels_N <- 100
curRun <- "runs12"
ROIdata <- focal_data_100_runs1234_VOEextract_bothdomains %>%
filter(ROI_name_final == curROI) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
}
else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID))
focal_region <- ROIdata$focal_region[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
focal_modelsummaries_100_VOEextract_bothdomains_init <-
rbind(focal_modelsummaries_100_VOEextract_bothdomains_init,
modelsummary)
}
focal_modelsummaries_100_VOEextract_bothdomains <- focal_modelsummaries_100_VOEextract_bothdomains_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
focal_modelsummaries_100_VOEextract_bothdomains$effect <-
as.factor(focal_modelsummaries_100_VOEextract_bothdomains$effect)
levels(focal_modelsummaries_100_VOEextract_bothdomains$effect) <-
c("intercept", "event")
rownames(focal_modelsummaries_100_VOEextract_bothdomains) <- NULL
# write_rds(focal_modelsummaries_100_VOEextract_bothdomains, here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))
focal_modelsummaries_100_VOEextract_bothdomains <- read_rds( here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))
DT::datatable(focal_modelsummaries_100_VOEextract_bothdomains, options = list(scrollX = TRUE))
options(contrasts = c("contr.sum", "contr.poly"))
single_betas0 <- readRDS(here("outputs/univariate_data/study2_univariate_data_singlebetas.Rds")) %>%
select(-ROI_category)
regioninfo_category <- regioninfo %>%
select(ROI_name_final, ROI_category)
single_betas <- left_join(single_betas0, regioninfo_category)
## Joining, by = "ROI_name_final"
single_betas$ROI_name_final <- as.factor(single_betas$ROI_name_final)
single_betas$event <- as.factor(single_betas$event)
single_betas$event <- relevel(single_betas$event, ref = "fam")
single_betas$task <- factor(single_betas$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))
single_betas_summary_domain <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_domain_singlebetas.Rds"))
single_betas_summary_task <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_task_singlebetas.Rds"))
single_betas_summary_scenario <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_scenario_singlebetas.Rds"))
single_betas_summary_scenario$event <- relevel(single_betas_summary_scenario$event, ref = "fam")
single_betas_focal <- single_betas %>%
filter(focal_region == 1)
# solidity - 25a88b
# permanence - 19626a
# goal - f24b00
# efficiency - deae1e
# agent-solidity - 6600cd
# missing-barrier - fba0d0
ggplot(
single_betas_focal %>%
filter(
str_length(extracted_run_number) == 4,
event != "fam"
),
aes(task, meanbeta, alpha = event, fill = task)
) +
geom_bar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
geom_errorbar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
scale_fill_manual(values = c(
"#25a88b",
"#19626a",
"#f24b00",
"#deae1e",
"#6600cd",
"#fba0d0"
)) + # ,
# geom_point() +
theme(axis.text.x = element_blank()) +
# scale_alpha_continuous(range = c(0.5, 1)) +
facet_wrap( ~ ROI_name_final, scales = "free", nrow = 2) +
ggtitle("Runs 1-4, Focal regions")
## Warning: Using alpha for a discrete variable is not advised.
ggplot(
single_betas_focal %>%
filter(ROI_name_final == "supramarginal_L"),
aes(
vis_statistics_video_group,
meanbeta,
alpha = event,
fill = task
)
) +
geom_bar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
geom_errorbar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
scale_fill_manual(values = c(
"#25a88b",
"#19626a",
"#f24b00",
"#deae1e",
"#6600cd",
"#fba0d0"
)) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
ggtitle("Region: supramarginal_L")
## Warning: Using alpha for a discrete variable is not advised.
Plotting each event type in each task against run number, per ROI.
ggplot(single_betas_focal, aes(as.integer(extracted_run_number), meanbeta, colour = event)) +
geom_smooth(method = "glm") +
facet_wrap( ~ ROI_name_final + task, nrow = 3, scales = "free") +
scale_colour_brewer(palette = "Set2")
## `geom_smooth()` using formula = 'y ~ x'
Visualizing visual features against responses.
ggplot(single_betas_focal,
aes(normalized_iqr_contrast, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_luminance, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_motion, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_hsf, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_lsf, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_rect, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_curv, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'
NOTE: Here we are using summed contrasts, which means that main effects are interpreted with respect to the grand mean.
options(contrasts = c("contr.sum", "contr.poly"))
single_betas_focal$event <- relevel(single_betas_focal$event, ref = "unexp")
single_betas_focal_psychphys_runs12 <- single_betas_focal %>%
filter(domain != "psychology-environment",
event != "fam") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")
ROIs <- unique(single_betas_focal_psychphys_runs12$ROI_name_final) %>% as.vector()
single_betas_focal_psychphys_runs12_summaries <- data.frame()
rownames_no_intercept <- c("event1", "domain1", "event:domain")
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group)
) # had to remove run random int (convergence))
} else {
ROI_category <- "general"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group)
)
}
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2), BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
ROI_focal <- ROIdata$focal_region[1]
cur_top_voxels_N <- 100
curRun <- "runs12"
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6,],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
single_betas_focal_psychphys_runs12_summaries <-
rbind(single_betas_focal_psychphys_runs12_summaries,
modelsummary)
}
single_betas_focal_psychphys_runs12_summaries <- single_betas_focal_psychphys_runs12_summaries %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_focal_psychphys_runs12_summaries$effect <-
as.factor(single_betas_focal_psychphys_runs12_summaries$effect)
levels(single_betas_focal_psychphys_runs12_summaries$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(single_betas_focal_psychphys_runs12_summaries) <- NULL
# write_rds(single_betas_focal_psychphys_runs12_summaries, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))
single_betas_focal_psychphys_runs12_summaries <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))
DT::datatable(dplyr::arrange(single_betas_focal_psychphys_runs12_summaries, effect), options = list(scrollX = TRUE))
single_betas_focal_visualstats_summary_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
} else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata, formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
ROI_focal <- ROIdata$focal_region[1]
cur_top_voxels_N <- 100
curRun <- "runs12"
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:12,],
effect_sizes
# BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized"
)
single_betas_focal_visualstats_summary_init <-
rbind(single_betas_focal_visualstats_summary_init,
modelsummary)
}
single_betas_focal_visualstats_summary <- single_betas_focal_visualstats_summary_init %>%
mutate_at(c(7:14), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_focal_visualstats_summary$effect <-
as.factor(single_betas_focal_visualstats_summary$effect)
levels(single_betas_focal_visualstats_summary$effect)[1:4] <-
c("intercept", "domain", "event", "event:domain")
rownames(single_betas_focal_visualstats_summary) <- NULL
# write_rds(single_betas_focal_visualstats_summary, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))
single_betas_focal_visualstats_summary <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))
DT::datatable(dplyr::arrange(single_betas_focal_visualstats_summary, effect), options = list(scrollX = TRUE))
V1_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "V1_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(V1_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "V1_bilateral")
##
## REML criterion at convergence: 3916
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.065 -0.581 -0.022 0.648 4.223
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 12.5313 3.540
## extracted_run_number (Intercept) 0.0802 0.283
## Residual 8.1251 2.850
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.0972 0.6683 25.6178 13.61 3.1e-13 ***
## event1 0.0655 0.1033 726.0307 0.63 0.5259
## domain1 0.4682 0.1767 726.0014 2.65 0.0082 **
## normalized_iqr_contrast -0.0388 0.1052 726.0004 -0.37 0.7119
## normalized_mean_luminance -1.2349 0.2240 726.0008 -5.51 4.9e-08 ***
## normalized_mean_motion 0.7183 0.1649 726.0038 4.36 1.5e-05 ***
## normalized_mean_hsf 0.3798 0.2439 726.0013 1.56 0.1198
## normalized_mean_rect -0.3466 0.2017 726.0002 -1.72 0.0861 .
## normalized_mean_curv 0.1841 0.1807 726.0021 1.02 0.3086
## event1:domain1 0.0613 0.1030 726.0002 0.60 0.5520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.003
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.062 -0.036 0.084
## nrmlzd_mn_l 0.048 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.040 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.005 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.016 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.042 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
check_collinearity(V1_model_singlebetas)
MT_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "MT_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(MT_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "MT_bilateral")
##
## REML criterion at convergence: 3098
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.229 -0.591 0.009 0.605 3.287
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 6.351 2.520
## extracted_run_number (Intercept) 0.016 0.126
## Residual 2.716 1.648
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.42973 0.45984 30.85320 11.81 5.6e-13 ***
## event1 0.04583 0.05970 726.04405 0.77 0.44295
## domain1 0.18927 0.10213 726.00118 1.85 0.06427 .
## normalized_iqr_contrast -0.53523 0.06080 725.99966 -8.80 < 2e-16 ***
## normalized_mean_luminance -0.25298 0.12949 726.00030 -1.95 0.05112 .
## normalized_mean_motion 0.35157 0.09535 726.00465 3.69 0.00024 ***
## normalized_mean_hsf -1.38930 0.14100 726.00108 -9.85 < 2e-16 ***
## normalized_mean_rect 0.23006 0.11660 725.99933 1.97 0.04887 *
## normalized_mean_curv -0.59928 0.10449 726.00211 -5.74 1.4e-08 ***
## event1:domain1 0.00941 0.05954 725.99932 0.16 0.87444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.002
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.052 -0.036 0.084
## nrmlzd_mn_l 0.040 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.034 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.004 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.013 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.035 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
check_collinearity(MT_model_singlebetas)
APC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "antParietal_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(APC_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "antParietal_bilateral")
##
## REML criterion at convergence: 3635
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.046 -0.535 0.015 0.520 4.165
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.3797 2.093
## extracted_run_number (Intercept) 0.0278 0.167
## Residual 5.7731 2.403
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.7034 0.4017 23.0478 9.22 3.4e-09 ***
## event1 0.2533 0.0870 726.0513 2.91 0.00372 **
## domain1 0.1349 0.1489 726.0023 0.91 0.36517
## normalized_iqr_contrast -0.2247 0.0887 726.0005 -2.54 0.01145 *
## normalized_mean_luminance -0.3695 0.1888 726.0013 -1.96 0.05073 .
## normalized_mean_motion 0.2846 0.1390 726.0063 2.05 0.04095 *
## normalized_mean_hsf -0.6989 0.2056 726.0022 -3.40 0.00071 ***
## normalized_mean_rect -0.0208 0.1700 726.0002 -0.12 0.90257
## normalized_mean_curv -0.4586 0.1523 726.0034 -3.01 0.00270 **
## event1:domain1 -0.0903 0.0868 726.0002 -1.04 0.29868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.004
## domain1 -0.005 -0.055
## nrmlzd_qr_c -0.087 -0.036 0.084
## nrmlzd_mn_l 0.067 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.056 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.007 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.022 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.058 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
check_collinearity(APC_model_singlebetas)
RFC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "precentral_inferiorfrontal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
summary(RFC_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "precentral_inferiorfrontal_R")
##
## REML criterion at convergence: 3629
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.369 -0.594 0.030 0.630 3.362
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 3.19 1.79
## Residual 5.81 2.41
## Number of obs: 768, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.8432 0.3324 32.8136 11.56 4.1e-13 ***
## event1 0.2549 0.0873 727.0000 2.92 0.0036 **
## domain1 0.2856 0.1493 727.0000 1.91 0.0562 .
## normalized_iqr_contrast -0.2156 0.0889 727.0000 -2.43 0.0155 *
## normalized_mean_luminance -0.2577 0.1893 727.0000 -1.36 0.1738
## normalized_mean_motion 0.1398 0.1394 727.0000 1.00 0.3161
## normalized_mean_hsf -0.6089 0.2062 727.0000 -2.95 0.0032 **
## normalized_mean_rect 0.1171 0.1705 727.0000 0.69 0.4926
## normalized_mean_curv -0.3032 0.1528 727.0000 -1.98 0.0475 *
## event1:domain1 -0.0143 0.0871 727.0000 -0.16 0.8696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.005
## domain1 -0.006 -0.055
## nrmlzd_qr_c -0.106 -0.036 0.084
## nrmlzd_mn_l 0.081 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.068 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.009 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.027 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.071 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
check_collinearity(RFC_model_singlebetas)
RSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(RSTS_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "superiortemporal_R")
##
## REML criterion at convergence: 3402
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.790 -0.573 0.014 0.589 3.773
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.50595 1.5830
## extracted_run_number (Intercept) 0.00173 0.0416
## Residual 4.29415 2.0722
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0826 0.2951 28.1989 7.06 1.1e-07 ***
## event1 0.1045 0.0751 726.1236 1.39 0.16439
## domain1 -0.3892 0.1284 726.0052 -3.03 0.00253 **
## normalized_iqr_contrast -0.1251 0.0765 726.0009 -1.64 0.10224
## normalized_mean_luminance -0.0367 0.1628 726.0027 -0.23 0.82180
## normalized_mean_motion -0.0225 0.1199 726.0149 -0.19 0.85101
## normalized_mean_hsf -0.6320 0.1773 726.0049 -3.56 0.00039 ***
## normalized_mean_rect 0.1491 0.1466 726.0000 1.02 0.30962
## normalized_mean_curv -0.0527 0.1314 726.0078 -0.40 0.68854
## event1:domain1 0.0298 0.0749 725.9999 0.40 0.69066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.005
## domain1 -0.005 -0.055
## nrmlzd_qr_c -0.103 -0.036 0.084
## nrmlzd_mn_l 0.078 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.066 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.009 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.026 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.069 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
check_collinearity(RSTS_model_singlebetas)
LSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + event_n_within_run +domain * event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(LSTS_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + event_n_within_run + domain * event +
## normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion +
## normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv +
## (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "superiortemporal_L")
##
## REML criterion at convergence: 3719
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.839 -0.511 0.058 0.546 3.911
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.386 2.094
## extracted_run_number (Intercept) 0.111 0.333
## Residual 6.414 2.533
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.1481 0.4877 11.9224 2.35 0.037 *
## event1 0.0625 0.0919 725.0360 0.68 0.496
## domain1 -0.3916 0.1570 724.9993 -2.50 0.013 *
## event_n_within_run -0.0305 0.0102 732.2838 -2.98 0.003 **
## normalized_iqr_contrast -0.0894 0.0934 724.9989 -0.96 0.339
## normalized_mean_luminance 0.1916 0.1990 725.0016 0.96 0.336
## normalized_mean_motion -0.0932 0.1465 725.0002 -0.64 0.525
## normalized_mean_hsf -0.3465 0.2167 725.0001 -1.60 0.110
## normalized_mean_rect -0.2323 0.1792 725.0003 -1.30 0.195
## normalized_mean_curv 0.0402 0.1607 725.0061 0.25 0.803
## event1:domain1 -0.0997 0.0916 725.0082 -1.09 0.276
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 evn___ nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m
## event1 -0.015
## domain1 0.000 -0.055
## evnt_n_wth_ -0.376 0.050 -0.010
## nrmlzd_qr_c -0.072 -0.037 0.084 -0.011
## nrmlzd_mn_l 0.066 -0.027 0.212 -0.021 -0.032
## nrmlzd_mn_m 0.047 0.086 -0.585 0.004 -0.381 -0.212
## nrmlzd_mn_h -0.001 -0.062 0.714 -0.014 0.399 0.344 -0.605
## nrmlzd_mn_r 0.013 -0.009 0.324 0.018 -0.017 -0.432 -0.089
## nrmlzd_mn_c -0.039 -0.056 0.442 -0.031 0.365 0.498 -0.615
## event1:dmn1 0.011 -0.006 0.030 -0.037 0.017 0.002 -0.049
## nrmlzd_mn_h nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## evnt_n_wth_
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r 0.296
## nrmlzd_mn_c 0.497 -0.426
## event1:dmn1 0.030 0.011 0.024
check_collinearity(LSTS_model_singlebetas)
RSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(RSMG_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "supramarginal_R")
##
## REML criterion at convergence: 3449
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.579 -0.562 0.005 0.565 6.309
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 7.3025 2.702
## extracted_run_number (Intercept) 0.0787 0.281
## Residual 4.3692 2.090
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.5859 0.5250 20.0189 4.93 0.000081 ***
## event1 0.1017 0.0757 726.0185 1.34 0.179
## domain1 0.2934 0.1295 726.0007 2.27 0.024 *
## normalized_iqr_contrast -0.1754 0.0771 726.0001 -2.27 0.023 *
## normalized_mean_luminance -0.4043 0.1642 726.0003 -2.46 0.014 *
## normalized_mean_motion 0.1410 0.1209 726.0021 1.17 0.244
## normalized_mean_hsf -0.4244 0.1788 726.0007 -2.37 0.018 *
## normalized_mean_rect 0.1570 0.1479 725.9999 1.06 0.289
## normalized_mean_curv -0.2804 0.1325 726.0011 -2.12 0.035 *
## event1:domain1 0.1810 0.0755 725.9999 2.40 0.017 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.003
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.058 -0.036 0.084
## nrmlzd_mn_l 0.044 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.037 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.005 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.015 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.039 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
plot(allEffects(RSMG_model_singlebetas))
lsmeans(RSMG_model_singlebetas, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.565 0.213 726 2.649 0.0080
##
## domain = psychology-action:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.158 0.214 726 -0.739 0.4600
##
## Degrees-of-freedom method: kenward-roger
check_collinearity(RSMG_model_singlebetas)
LSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(LSMG_model_singlebetas)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "supramarginal_L")
##
## REML criterion at convergence: 3379
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.060 -0.600 0.012 0.579 5.640
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.11543 2.0287
## extracted_run_number (Intercept) 0.00379 0.0616
## Residual 4.07805 2.0194
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.4037 0.3715 30.2024 6.47 3.7e-07 ***
## event1 0.0665 0.0731 726.1058 0.91 0.3636
## domain1 0.2683 0.1252 726.0045 2.14 0.0324 *
## normalized_iqr_contrast -0.0619 0.0745 726.0008 -0.83 0.4064
## normalized_mean_luminance -0.4990 0.1587 726.0024 -3.14 0.0017 **
## normalized_mean_motion -0.0931 0.1168 726.0128 -0.80 0.4256
## normalized_mean_hsf -0.0259 0.1728 726.0042 -0.15 0.8808
## normalized_mean_rect 0.2933 0.1429 726.0001 2.05 0.0405 *
## normalized_mean_curv -0.0639 0.1280 726.0067 -0.50 0.6176
## event1:domain1 0.2042 0.0730 726.0000 2.80 0.0053 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.004
## domain1 -0.004 -0.055
## nrmlzd_qr_c -0.079 -0.036 0.084
## nrmlzd_mn_l 0.061 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.051 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.007 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.020 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.053 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023
plot(allEffects(LSMG_model_singlebetas))
lsmeans(LSMG_model_singlebetas, pairwise ~ event | domain)$contrasts
## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.541 0.206 726 2.626 0.0088
##
## domain = psychology-action:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.275 0.207 726 -1.330 0.1838
##
## Degrees-of-freedom method: kenward-roger
check_collinearity(LSMG_model_singlebetas)
single_betas_pertask_data <- single_betas %>%
filter(event != "fam") %>%
filter(str_length(extracted_run_number) == 4)
tasks <- unique(single_betas_pertask_data$task)
ROIs <- unique(single_betas_pertask_data$ROI_name_final)
single_betas_pertask_data$event <- relevel(single_betas_pertask_data$event, ref = "unexp")
runs <- c("all_runs", "runs12")
single_betas_pertask_results <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (task in tasks) {
cur_task <- task
for (run in runs) {
if (run == "all_runs") {
ROIdata <-
single_betas_pertask_data %>% filter(ROI_name_final == curROI,
task == cur_task)
}
else if (run == "runs12") {
ROIdata <-
single_betas_pertask_data %>%
filter(ROI_name_final == curROI,
task == cur_task) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
}
# print("Current ROI: ", curROI)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID)) # had to remove run random int (convergence))
focal_region <- ROIdata$focal_region[1]
task_domain <- as.character(ROIdata$domain[1])
ROI_category <- ROIdata$ROI_category[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes <- standardize_parameters(model) %>%
as.data.frame() %>%
select(Std_Coefficient)
cur_top_voxels_N <- 100
curRun <- run
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
cur_task,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"task",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"B_standardized",
"BF",
"BF_interpretation"
)
single_betas_pertask_results <-
rbind(single_betas_pertask_results, modelsummary)
}
}
}
single_betas_pertask_results <- single_betas_pertask_results %>%
mutate_at(c(8:16), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_pertask_results$effect <-
as.factor(single_betas_pertask_results$effect)
levels(single_betas_pertask_results$effect) <-
c("intercept", "event")
rownames(single_betas_pertask_results) <- NULL
# write_rds(single_betas_pertask_results, here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
DT::datatable(dplyr::arrange(single_betas_pertask_results, extracted_run_number, effect), options = list(scrollX = TRUE))
overlap_data_prelim <- read_csv(here("input_data/ROI_overlap_outputs_Apr19_2023.csv")) %>%
rename(ROI_2 = roi_2,
subjID = participantID)
## Parsed with column specification:
## cols(
## participantID = col_character(),
## ROI_1 = col_character(),
## roi_2 = col_character(),
## overlap_dices_coef = col_double()
## )
MD_physics_APC_bilateral <- overlap_data_prelim %>%
filter((
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
ROI_2 == "MD_ROI_3_13_bilateral"
) |
(
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
ROI_2 == "MD_ROI_3_13_bilateral"
)
) %>%
mutate(hemisphere = "bilateral")
MD_physics_APC_byhem <- overlap_data_prelim %>%
filter((
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
ROI_2 == "MD_ROI_13"
) |
(
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
ROI_2 == "MD_ROI_3"
)
) %>%
mutate(hemisphere = case_when(ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" ~ "right",
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" ~ "left"))
names(MD_physics_APC_bilateral)
## [1] "subjID" "ROI_1" "ROI_2"
## [4] "overlap_dices_coef" "hemisphere"
names(MD_physics_APC_byhem)
## [1] "subjID" "ROI_1" "ROI_2"
## [4] "overlap_dices_coef" "hemisphere"
MD_physics_APC_bilateral <- rbind(MD_physics_APC_bilateral, MD_physics_APC_byhem)
overlap_summary <- ggplot(MD_physics_APC_bilateral, aes(hemisphere, overlap_dices_coef)) +
geom_boxplot() +
geom_point(fun = "mean", stat = "summary", shape = 18, colour = "red", size = 4) +
geom_jitter(height = 0, width = .1) +
ylab("Dice's Coefficient") +
xlab("APC ROI (Bilateral = Confirmatory)")
# geom_text_repel(aes(label = subjID), size = 3) +
# ggtitle("Dice's Coefficient between LSMG and RSMG \n to bilateral APC, left APC, and right APC")
overlap_summary
APC_SMG_overlap_results <- MD_physics_APC_bilateral %>%
group_by(hemisphere) %>%
summarise(mean_DC = mean(overlap_dices_coef),
median_DC = median(overlap_dices_coef),
range_low = range(overlap_dices_coef)[1],
range_high = range(overlap_dices_coef)[2])
MD_physics_APC_bilateral %>%
filter(overlap_dices_coef == 0) %>%
group_by(hemisphere) %>%
summarise(n = n())
MD_physics_APC_bilateral %>%
filter(overlap_dices_coef != 0,
hemisphere == "bilateral") %>%
group_by(subjID)
MD_physics_APC_byhem %>%
group_by(hemisphere) %>%
summarise(mean_DC = mean(overlap_dices_coef),
median_DC = median(overlap_dices_coef),
range_low = range(overlap_dices_coef)[1],
range_high = range(overlap_dices_coef)[2])
MD_physics_APC_byhem %>%
filter(overlap_dices_coef == 0) %>%
group_by(hemisphere) %>%
summarise(n = n())
# write_rds(APC_SMG_overlap_results, path = here("outputs/univariate_results/APC_SMG_overlap_results.Rds"))
single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
# filter((domain == "both" & extracted_run_number == "all_runs") | (domain != "both" & extracted_run_number == "runs12"))
tasks <- c("infer-constraint", "agent-solidity")
single_betas_pertask_results_focal <- single_betas_pertask_results %>%
filter(focal_region == 1,
effect == "event") %>%
filter((extracted_run_number == "all_runs" & task %in% tasks) | (extracted_run_number == "runs12" & !task %in% tasks)) %>%
select(extracted_run_number, ROI, task, domain, B_standardized, BF, B, SE, CI_95_lower, CI_95_upper, star)
single_betas_pertask_results_focal$task <- factor(single_betas_pertask_results_focal$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))
single_betas_pertask_results_focal$ROI <- factor(single_betas_pertask_results_focal$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral"))
pertask_b_pertask <- ggplot(single_betas_pertask_results_focal,
aes(ROI, B, colour = task, group = task, label = star)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
geom_text(size = 4, aes(ROI, B + SE + .02)) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ task, nrow = 1) +
scale_colour_manual(values = c("#25a88b", "#19626a", "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.6, .8))
pertask_b_perROI <- ggplot(single_betas_pertask_results_focal,
aes(task, B, colour = task, group = task, label = star)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
geom_text(size = 4, aes(task, B + SE + .02)) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ ROI, nrow = 1) +
scale_colour_manual(values = c("#25a88b", "#19626a", "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.6, .8))
pertask_b_pertask
pertask_b_perROI
This analysis tests for the tolerance of the primary VOE effects to low level visual confounds.
NOTE: Here we are using treatment contrasts, which means that main effects are interpreted with respect to the reference level.
options(contrasts = c("contr.treatment", "contr.poly"))
data_for_modeling <- single_betas %>%
filter(event != "fam",
top_n_voxels == "100",
domain != "psychology-environment") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2") %>%
filter(!str_detect(ROI_name_final, "bilateral"))
data_for_modeling$ROI_name_final <- as.factor(data_for_modeling$ROI_name_final)
ROIs <- sort(unique(data_for_modeling$ROI_name_final))
data_for_modeling$event <- relevel(data_for_modeling$event, ref = "exp")
data_for_modeling$domain <- relevel(data_for_modeling$domain, ref = "psychology-action")
modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
# FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
physics_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "physics")
psychology_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "psychology-action")
# had to remove random intercept for runs (convergence))
physics_model <- lmer(data = physics_data,
formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
psychology_model <- lmer(data = psychology_data,
formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
physics_event_effect_size <-
standardize_parameters(physics_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
psychology_event_effect_size <- standardize_parameters(psychology_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
# SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
expected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "exp")
unexpected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "unexp")
# had to remove random intercept for runs (convergence))
exp_model <- lmer(data = expected_data,
formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
unexp_model <- lmer(data = unexpected_data,
formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
exp_domain_effect_size <-
standardize_parameters(exp_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
unexp_domain_effect_size <-
standardize_parameters(unexp_model) %>%
as.data.frame() %>%
select(Std_Coefficient) %>%
slice(2)
if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
ROI_domain <- "general"
} else {
ROI_domain <- "specific"
}
modelsummary <-
cbind(
ROI_domain,
physics_data$ROI_category[1],
curROI,
physics_event_effect_size,
psychology_event_effect_size,
exp_domain_effect_size,
unexp_domain_effect_size
) %>%
as.data.frame()
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}
colnames(modelsummaries_init) <-
c("ROI_domain",
"ROI_category",
"ROI",
"physics_event_effect_size",
"psychology_event_effect_size",
"exp_domain_effect_size",
"unexp_domain_effect_size"
)
rownames(modelsummaries_init) <- NULL
# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))
univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))
plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size, physics_event_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-0.2, .75), ylim = c(-0.2, .75)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
xlab("VOE Effect Size (d), Psychology Events") +
ylab("VOE Effect Size (d), Physics Events") +
ggtitle(label = "Exp 2 VOE effects for each domain across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
xlab("Domain Effect Size (d), Expected Events") +
ylab("Domain Effect Size (d), Unexpected Events") +
ggtitle(label = "Exp 2 Domain effects for each event across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot1
## `geom_smooth()` using formula = 'y ~ x'
plot2
## `geom_smooth()` using formula = 'y ~ x'
# instead of testing whether the linear relationship between x and y is 0,
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
group_by(ROI_domain) %>%
summarise(
cor_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value,
cor_event_within_domain =
np.cor.test(
psychology_event_effect_size,
physics_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_event_within_domain =
np.cor.test(
psychology_event_effect_size,
physics_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value
)
observed_cors_ind
# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study2_univar_manyregions_observed_cor.Rds"))
set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
data <- data[indices,]
cor1 <-
np.cor.test(DS_results$exp_domain_effect_size,
DS_results$unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
cor2 <-
np.cor.test(
data$psychology_event_effect_size,
data$physics_event_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
return(cor1 - cor2)
}
DS_results <- univar_manyregions_results %>%
filter(ROI_domain == "specific")
doMC::registerDoMC(cores = 2) # for running in parallel
# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
# data = DS_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i"
# )
# saveRDS(res_boot_DS, here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))
res_boot_DS <- read_rds(here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))
## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))
# saveRDS(ds_results, here("outputs/univariate_results/study2_univariate_dsregions_summary.Rds"))
DG_results <- univar_manyregions_results %>%
filter(ROI_domain == "general")
## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i")
# saveRDS(res_boot_DG, here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))
plot(res_boot_DG)
## Retrieve the empirical 95% confidence interval (adjusted to 90% because of one-tailed prediction):
# dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study2_univariate_dgregions_summary.Rds"))
V1_bilat_famvstest <- focal_data %>%
filter(ROI_name_final == "V1_bilateral",
domain != "both",
top_n_voxels == 100) %>%
filter(str_length(extracted_run_number) == 4) %>% # restrict analysis to same as confirm
mutate(event_famtest = case_when(event == "fam" ~ "fam",
TRUE ~ "test")) # compare fam to test
MT_bilat_famvstest <- focal_data %>%
filter(ROI_name_final == "MT_bilateral",
domain != "both",
top_n_voxels == 100) %>%
filter(str_length(extracted_run_number) == 4) %>% # restrict analysis to same as confirm
mutate(event_famtest = case_when(event == "fam" ~ "fam",
TRUE ~ "test")) # compare fam to test
nonvis_famvstest <- focal_data %>%
filter(ROI_category != "early_visual",
domain != "both",
top_n_voxels == 100) %>%
filter(str_length(extracted_run_number) == 4) %>% # restrict analysis to same as confirm analysis to same as confirm
mutate(event_famtest = case_when(event == "fam" ~ "fam",
TRUE ~ "test")) # compare fam to test
focalregions_famvstest <- rbind(V1_bilat_famvstest, MT_bilat_famvstest, nonvis_famvstest)
ggplot(focalregions_famvstest,
aes(event_famtest, meanbeta, fill = ROI_category)) +
geom_boxplot() +
facet_wrap(~ ROI_category,
nrow = 1,
scales = "free") +
scale_fill_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
stat_summary(geom = "errorbar", fun.data = "mean_cl_boot", colour = "white", width = .2) +
ylab("Average amplitude of response") +
xlab("Event Type")
V1_visnovelty_model <- lmer(data = V1_bilat_famvstest, formula = meanbeta ~ event_famtest + (1 | subjectID))
summary(V1_visnovelty_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event_famtest + (1 | subjectID)
## Data: V1_bilat_famvstest
##
## REML criterion at convergence: 3530
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.907 -0.590 0.029 0.623 3.239
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 12.77 3.57
## Residual 5.07 2.25
## Number of obs: 762, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.660 0.647 33.076 14.92 3e-16 ***
## event_famtesttest -0.351 0.173 729.009 -2.03 0.043 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## evnt_fmtstt -0.178
plot(allEffects(V1_visnovelty_model))
check_model(V1_visnovelty_model)
MT_visnovelty_model <- lmer(data = MT_bilat_famvstest, formula = meanbeta ~ event_famtest + (1 | subjectID))
summary(MT_visnovelty_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event_famtest + (1 | subjectID)
## Data: MT_bilat_famvstest
##
## REML criterion at convergence: 2662
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.745 -0.593 0.030 0.630 3.581
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 6.63 2.57
## Residual 1.58 1.26
## Number of obs: 762, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.3730 0.4619 32.2508 11.63 4.4e-13 ***
## event_famtesttest -0.2770 0.0967 729.0060 -2.86 0.0043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## evnt_fmtstt -0.140
check_model(V1_visnovelty_model)
V1_visnovelty_model_results <- cbind(gen.m(V1_visnovelty_model), gen.ci(V1_visnovelty_model)[3:4,])
## Computing profile confidence intervals ...
MT_visnovelty_model_results <- cbind(gen.m(MT_visnovelty_model), gen.ci(MT_visnovelty_model)[3:4,])
## Computing profile confidence intervals ...
knitr::kable(V1_visnovelty_model_results)
| B | SE | df | t | p | CI_95_lower | CI_95_upper | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 9.660 | 0.647 | 33.1 | 14.92 | 0.000 | 8.374 | 10.945 |
| event_famtesttest | -0.351 | 0.173 | 729.0 | -2.03 | 0.043 | -0.691 | -0.012 |
knitr::kable(MT_visnovelty_model_results)
| B | SE | df | t | p | CI_95_lower | CI_95_upper | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 5.373 | 0.462 | 32.3 | 11.63 | 0.000 | 4.455 | 6.291 |
| event_famtesttest | -0.277 | 0.097 | 729.0 | -2.86 | 0.004 | -0.467 | -0.087 |
# write_rds(V1_visnovelty_model_results, path = here("outputs/univariate_results/V1_visnovelty_results.Rds"))
# write_rds(MT_visnovelty_model_results, path = here("outputs/univariate_results/MT_visnovelty_results.Rds"))
allregions_visnovelty_model <- lmer(data = focalregions_famvstest,
formula = meanbeta ~ ROI_category * event_famtest + (1 | subjectID) + (1|ROI_name_final))
summary(allregions_visnovelty_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ ROI_category * event_famtest + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: focalregions_famvstest
##
## REML criterion at convergence: 28513
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.014 -0.579 -0.047 0.509 5.253
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.48 1.22
## ROI_name_final (Intercept) 2.53 1.59
## Residual 6.13 2.48
## Number of obs: 6096, groups: subjectID, 32; ROI_name_final, 8
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 7.511 1.150 4.349 6.53
## ROI_categoryMD -5.023 1.598 4.051 -3.14
## ROI_categoryphysics -6.086 1.598 4.051 -3.81
## ROI_categorypsychology -6.970 1.598 4.051 -4.36
## event_famtesttest -0.314 0.135 6053.003 -2.33
## ROI_categoryMD:event_famtesttest 1.167 0.190 6053.003 6.13
## ROI_categoryphysics:event_famtesttest 1.043 0.190 6053.003 5.48
## ROI_categorypsychology:event_famtesttest 0.904 0.190 6053.003 4.75
## Pr(>|t|)
## (Intercept) 0.0021 **
## ROI_categoryMD 0.0342 *
## ROI_categoryphysics 0.0185 *
## ROI_categorypsychology 0.0117 *
## event_famtesttest 0.0196 *
## ROI_categoryMD:event_famtesttest 9.1e-10 ***
## ROI_categoryphysics:event_famtesttest 4.4e-08 ***
## ROI_categorypsychology:event_famtesttest 2.1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ROI_cMD ROI_ctgryph ROI_ctgryps evnt_f ROI_MD:
## ROI_ctgryMD -0.695
## ROI_ctgryph -0.695 0.500
## ROI_ctgryps -0.695 0.500 0.500
## evnt_fmtstt -0.078 0.056 0.056 0.056
## ROI_ctgMD:_ 0.055 -0.079 -0.040 -0.040 -0.707
## ROI_ctgryph:_ 0.055 -0.040 -0.079 -0.040 -0.707 0.500
## ROI_ctgryps:_ 0.055 -0.040 -0.040 -0.079 -0.707 0.500
## ROI_ctgryph:_
## ROI_ctgryMD
## ROI_ctgryph
## ROI_ctgryps
## evnt_fmtstt
## ROI_ctgMD:_
## ROI_ctgryph:_
## ROI_ctgryps:_ 0.500
plot(allEffects(allregions_visnovelty_model))
allregions_visnovelty_model_simpler <- lmer(data = focalregions_famvstest,
formula = meanbeta ~ ROI_category + event_famtest + (1 | subjectID) + (1|ROI_name_final))
anova(allregions_visnovelty_model_simpler, allregions_visnovelty_model)
## refitting model(s) with ML (instead of REML)
results <- lsmeans(allregions_visnovelty_model, pairwise ~ event_famtest | ROI_category, pbkrtest.limit = 6096)$contrasts %>% as.data.frame()
results
allregions_visnovelty_model_results <- cbind(gen.m(allregions_visnovelty_model), gen.ci(allregions_visnovelty_model)[3:10,])
## Computing profile confidence intervals ...
knitr::kable(allregions_visnovelty_model_results)
| B | SE | df | t | p | CI_95_lower | CI_95_upper | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 7.511 | 1.150 | 4.35 | 6.53 | 0.002 | 2.432 | 2.52 |
| ROI_categoryMD | -5.023 | 1.598 | 4.05 | -3.14 | 0.034 | 5.670 | 9.35 |
| ROI_categoryphysics | -6.086 | 1.598 | 4.05 | -3.81 | 0.019 | -7.571 | -2.48 |
| ROI_categorypsychology | -6.970 | 1.598 | 4.05 | -4.36 | 0.012 | -8.634 | -3.54 |
| event_famtesttest | -0.314 | 0.135 | 6053.00 | -2.33 | 0.020 | -9.518 | -4.42 |
| ROI_categoryMD:event_famtesttest | 1.167 | 0.190 | 6053.00 | 6.13 | 0.000 | -0.578 | -0.05 |
| ROI_categoryphysics:event_famtesttest | 1.043 | 0.190 | 6053.00 | 5.48 | 0.000 | 0.794 | 1.54 |
| ROI_categorypsychology:event_famtesttest | 0.904 | 0.190 | 6053.00 | 4.75 | 0.000 | 0.671 | 1.42 |
# write_rds(allregions_visnovelty_model_results, path = here("outputs/univariate_results/allfocalregions_visnovelty_modelsummary.Rds"))
# write_rds(results, path = here("outputs/univariate_results/allfocalregions_visnovelty_results_per_ROI_cat.Rds"))